function results = deterministic_mimo_capacity(varargin)
% 确定性MIMO信道容量分析
% 输入参数:
%   'SNR_dB' - SNR范围 (dB), 默认: -5:5:30
%   'H_det' - 确定性信道矩阵, 默认: 4x4随机信道
% 输出:
%   results - 包含容量结果、奇异值、功率分配等信息的结构体

% 解析输入参数
p = inputParser;
addParameter(p, 'SNR_dB', -5:5:30);
addParameter(p, 'H_det', sqrt(0.5) * (randn(4, 4) + 1i * randn(4, 4)));
parse(p, varargin{:});

SNR_dB = p.Results.SNR_dB;
H_det = p.Results.H_det;
SNR_linear = 10.^(SNR_dB/10);

% 添加路径
addpath('../Common');

% 获取颜色定义
colors = color_definitions();

fprintf('=== 确定性MIMO信道容量分析 ===\n');

% 信道矩阵分析
fprintf('确定性信道矩阵 (4x4):\n');
fprintf('Frobenius范数: %.2f\n', norm(H_det, 'fro'));
fprintf('条件数: %.2f\n', cond(H_det));

% 计算不同SNR下的容量
det_capacity_csi = zeros(length(SNR_dB), 1);
det_capacity_no_csi = zeros(length(SNR_dB), 1);

for snr_idx = 1:length(SNR_dB)
    snr = SNR_linear(snr_idx);
    
    % 发送端已知CSI (注水算法)
    [U, S, V] = svd(H_det);
    singular_values = diag(S);
    det_capacity_csi(snr_idx) = water_filling_capacity(singular_values, snr);
    
    % 发送端未知CSI (等功率分配)
    det_capacity_no_csi(snr_idx) = log2(det(eye(4) + (snr/4) * H_det * H_det'));
end

% 绘制确定性信道容量
figure('Name', '确定性MIMO信道容量', 'Position', [100, 100, 1200, 800]);

subplot(2,2,1);
plot(SNR_dB, det_capacity_csi, 'b-', 'LineWidth', 2);
hold on;
plot(SNR_dB, det_capacity_no_csi, 'r--', 'LineWidth', 2);
grid on;
xlabel('SNR (dB)');
ylabel('容量 (bps/Hz)');
title('确定性4x4 MIMO信道容量');
legend('已知CSI', '未知CSI');

% 奇异值分析
[U, S, V] = svd(H_det);
singular_values = diag(S);

subplot(2,2,2);
plot(singular_values, 'o-', 'LineWidth', 2);
grid on;
xlabel('模式索引');
ylabel('奇异值');
title('信道奇异值分布');

% 功率分配 (特定SNR)
snr_fixed = 15; % dB
snr_linear = 10^(snr_fixed/10);
[opt_powers, water_level] = water_filling_powers(singular_values, snr_linear);

subplot(2,2,3);
bar(opt_powers);
grid on;
title(sprintf('注水功率分配 (SNR=%d dB)', snr_fixed));
xlabel('模式索引');
ylabel('分配功率');

% 容量贡献分析
mode_contributions = zeros(length(singular_values), 1);
for i = 1:length(singular_values)
    if opt_powers(i) > 0
        mode_contributions(i) = log2(1 + snr_linear * opt_powers(i) * singular_values(i)^2);
    end
end

subplot(2,2,4);
pie(mode_contributions);
title('各模式容量贡献');
legend(arrayfun(@(x) sprintf('模式%d', x), 1:length(singular_values), 'UniformOutput', false));

% 组织结果
results.SNR_dB = SNR_dB;
results.capacity_csi = det_capacity_csi;
results.capacity_no_csi = det_capacity_no_csi;
results.singular_values = singular_values;
results.opt_powers = opt_powers;
results.water_level = water_level;
results.mode_contributions = mode_contributions;

end

function capacity = water_filling_capacity(singular_values, snr)
    % 注水算法容量计算
    num_modes = length(singular_values);
    total_power = length(singular_values); % 归一化总功率
    noise_power = 1/snr;
    
    % 注水功率分配
    water_level = water_filling_mimo(singular_values.^2, total_power, noise_power);
    
    % 计算容量
    capacity = 0;
    for i = 1:num_modes
        power_alloc = max(0, water_level - noise_power/singular_values(i)^2);
        if power_alloc > 0
            capacity = capacity + log2(1 + snr * power_alloc * singular_values(i)^2);
        end
    end
end

function [opt_powers, water_level] = water_filling_powers(singular_values, snr)
    % 注水功率分配
    num_modes = length(singular_values);
    total_power = num_modes;
    noise_power = 1/snr;
    
    % 注水算法
    water_level = water_filling_mimo(singular_values.^2, total_power, noise_power);
    
    % 功率分配
    opt_powers = zeros(num_modes, 1);
    for i = 1:num_modes
        opt_powers(i) = max(0, water_level - noise_power/singular_values(i)^2);
    end
end

function water_level = water_filling_mimo(eigenvalues, total_power, noise_power)
    % MIMO注水算法
    num_modes = length(eigenvalues);
    
    % 排序特征值
    [sorted_eig, sort_idx] = sort(eigenvalues, 'descend');
    
    % 寻找最优水位线
    water_level = 0;
    
    for k = 1:num_modes
        temp_level = (total_power + noise_power * sum(1./sorted_eig(1:k))) / k;
        temp_powers = max(0, temp_level - noise_power./sorted_eig(1:k));
        
        if all(temp_powers >= 0)
            water_level = temp_level;
        else
            break;
        end
    end
end